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o : 

The use of excessively long timesteps in dissipative particle dynamics simulations may produce 
simulation artifacts due to the generation of configurations which are not representative of the desired 
| canonical ensemble. The configurational temperature, amongst other quantities, may be used to 

. assess the extent of the deviation from equilibrium. This paper presents results for simulations of 

' models of water, and lipid bilayer membranes, to illustrate the nature of the problems. 

I. INTRODUCTION 

Dissipative particle dynamics 0,0 (DPD) has become a popular tool for simulating the behaviour of both simple 
. and complex fluids. In outline, it consists of the solution of the classical equations of motion for a system of interacting 
particles, together with a set of stochastic and dissipative forces which control the temperature and allow one to choose 
the viscosity. It is usual to employ a very simple, repulsive, pair potential, often choosing the force law derived from it 
to vary linearly with separation within a specified cutoff range. As well as these conservative forces, the stochastic and 
dissipative forces also act in a pairwise fashion, so as to conserve momentum and ensure hydrodynamic behaviour. 
The particles represent fluid regions, rather than individual atoms and molecules: the softness and simplicity of 
the interactions permit the use of a long time step, compared with conventional molecular dynamics. This, and the 
acceleration of physical processes compared with those seen in more realistic simulations, gives an advantage of several 
orders of magnitude, at the cost of a very rough mapping onto specific molecular properties. The way to relate the 
)_pa 
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By tuning the interactions to represent hydrophilic and hydrophobic units, and linking together the particles in 
a simple bead-spring fashion, the DPD method has been applied to lipid bilayer membranes 4]. Self-assembly and 
phase diagrams have been investigated |£g, stress profiles determined p, and the effects of small molecules Q and 
t-H , membrane-bound proteins |8| studied. 

Several pap ers have attempted to optimise the computational algorithm used to solve the DPD equations 0, 0, 
ITU H^ . IT^L Il4l ITEL lill ITR Il8l | . The early emphasis was on the correct handling of the velocity-dependent dissipative 
forces, and ensuring the correct balance of these with the random forces, at an appropriate timestep, so as to generate 
the canonical ensemble |19| . Accordingly, most attention Q has focused on maintaining a kinetic temperature Tk 
(defined in terms of the kinetic energy per degree of freedom, via the equipartition theorem) reasonably close to the 
desired temperature T. A detailed analysis of the ideal fluid limit |20j. as a function of timestep, has been carried 
out. In addition, it has become standard practice to check the pair distribution function and diffusion coefficient for 
dependence on the time step, as diagnostic tests. 

Recently, increased attention has been paid to the suitability or otherwise of the time step for the integration of 
the conservative forces. This is particularly important in simulations of multi-component systems, polymers, and 
membranes, where some of the interaction potentials may be stronger than others, including intramolecular bond- 
bending and stretching terms. The usual criterion of molecular dynamics, energy conservation in the absence of 
random and dissipative forces, is not usually employed, and any problems may be disguised by the thermostatting. 
O 1 In order to generate meaningful results, it is essential to sample configurations from the canonical (Boltzmann) 
distribution, and any indicator of significant deviations must be taken seriously. Den Otter and Clarke pH were 
the first to point out that DPD, at the usual timesteps and state points, generates configurations for which the 
configurational temper ature T c differs significantly from the desired temperature T and the kinetic temperature T^. 
iDen Otter and Clarkd estimated this quantity by calculating (for small timesteps) the average potential energy (U) 
as a function of temperature in the vicinity of T, and fitting the results to a quadratic curve; then, this curve was 
used to convert averages of (U) measured at larger timesteps into an estimate of T c . Here a different, and more 
straightforward definition is used. The configurational temperature T c is written 

kBT < = ^Wr ■ (1) 

where fee is Boltzmann's constant, Vj the gradient, and V| the laplacian, of the potential energy U, with respect to 
the position of particle j. The summations may be over all particles in the system, or restricted to a single species, 
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or indeed to an individual particle. The expression has long been known as a hypervirial relation |22j,|23j. It was first 
introduced by Rugh [24[ as an independent estimate of temperature in simulations, and recommended as a diagnostic 
test for lack of equilibrium by Butler et al. |2fij . 

In the context of membrane simulations, Jakobsen et al. [26j observed significant differences in the kinetic tempera- 
tures of different species in DPD simulations, and recommended this as a simple diagnostic test of incorrect integration 
of the equations of motion. They also mentioned the possibility of monitoring the configurational temperature, but 
did not report any results for this. However, they observed several other simulation artifacts, including time-step de- 
pendence of the pressure profiles and particle densities within the membrane. Hafskjold et al. 27] have reemphasized 
that the timesteps typically recommended and used in DPD are too long to properly handle the discontinuity at the 
cutoff of the usual conservative force law: they investigated the possibility of softening this cutoff via spline functions, 
but no dramatic improvements at long time step were seen, and they recommended the use of shorter time steps. 

This paper reports simulation averages of the configurational temperature for a standard single-component DPD 
fluid, over a range of timesteps, using a variety of proposed DPD integration methods. The results show that, although 
some of these methods dramatically improve the control of the kinetic temperature Tk, the same improvements are 
not seen for the configurational temperature: deviations of T c from the desired values remain significant in all cases, 
at the typically recommended time steps. Results are then reported for the same membrane model used in Ref. |26| : 
largely as a result of intramolecular bonding terms in the potential, the deviations in T r are far worse than in the 
simple fluid case, and may be correlated with some of the other artifacts observed before |26j . The general conclusion 
is that the configurational temperature should be added to the list of diagnostic tests applied to DPD simulations, 
and that shorter timesteps should be employed in membrane simulations than hitherto. 



II. DYNAMICAL ALGORITHMS 

The DPD equations of motion for a simple fluid may be written P, 0, 0] 

dn = (pi/m)dt (2a) 

d Pi = E f% dt + tf dt + d pl ■ ( 2b ) 

The particles are all assumed to have the same mass m. The conservative forces usually take the form 

fij=^(rij)fij (3a) 

I 1 — t/Tc r < r c ,„, . 

with u> (r) = { 1 ~ ■ ( 3b ) 
10 r > r c 

Here = — rj, rij — r,j , r,j — Vij/rij. The parameter a determines the strength of the conservative interactions, 
and r c is the cutoff. The dissipative forces and random impulses are written 

fPj = -l^{n 3 f{v tj -fi^rij 
dp* = <TU}(n j )<r ij dW ij , 

where Vij — Vi — Vj, and pi = mVi. A choice has been made here to use the same weighting function u)(r) in 
the specification of conservative, dissipative, and random terms. The dissipative friction 7 is related to the random 
impulse strength a by the fluctuation-dissipation theorem a 2 = . The quantity dWij is the time derivative of 

a Wiener process: over a timestep At, AWij is chosen from a normal distribution with zero mean and variance At, 
and AWji = — AWij. Units are chosen such that the particle mass m, cutoff r c , and temperature T are all unity. 
Four algorithms will be discussed in this paper. The "standard" DPD algorithm consists of the following steps 
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Pi ■= Pi + /§At + /£ At + Ap* (4a) 



f*:=ft + j£/<iAt + /£At + Apg (4b) 
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:= n + (pi/m)At (4c) 
Calculate /§({rj), /°({rj, {&}), Apg({r<}) (4d) 

Calculate /° ({n}, (4f) 

Note the auxiliary momenta pi, used in the calculation of dissipative forces, and predicted using a parameter A. 
The value A = 0.5 conveniently makes pi — pi, but A = 0.65 was found by Groot and Warren [3| to be optimal 
for the particular fluid that will be studied here. This paper will use A = 0.65. Note also the second computation 
of dissipative forces, step (UJJ, introduced by Gibson and Chen Iteration of the last two steps would lead to a 
self-consistent algorithm [12| not studied here. 

The second algorithm is the so-called "Scheme II" of Peters At each time step, the positions and momenta 
are updated using the conservative forces only: 

Vi :=n+ (pi/m)At (5b) 
Calculate /£({rj) (5c) 

Pi--=Pi+lJ2fiJ At - ( 5d ) 

Following this, for every pair selected in random order, the following update is performed: 

Pi :=Pi + /£At + Apg 
Pi ■= Pj - f$At - Ap]] 
with /° = -Oij(vij ■ fij)fij 
and ApJ* = 6,,Ail , ; /-,, . 

The parameters a,j and 6y are determined by T and 7 (or, equivalently, T and a) and by the timestep. As explained 
by Peters 0|j this generates the Maxwell-Boltzmann distribution of velocities and, at short time step, solves the 
original DPD equations. 

The third algorithm is due to Lowe 01 • This is not a solution of the original DPD equations, but rather a pairwise 
thermostat applied to a conventional integration algorithm for the conservative forces. As for the Peters method, at 
each step the positions and momenta are updated using steps l(5^|) - (|5dp above, using the conservative forces only. 
Then, examining every pair (in random order), with probability P — TAt the momenta are updated as follows: 

Pi :=Pi + Apfj 
Pi -=Pj- A Pfj 
with Apfj = -m 



C^/2k B T/m - (vij ■ fij) 



where £ is selected from a Gaussian distribution with zero mean and unit variance. This procedure reselects the relative 
velocity from the Maxwell-Boltzmann distribution. The key parameter is the stochastic randomization frequency P. 

The fourth algorithm to be considered here has been proposed by Stoyanov and Groot 0] . It is identical with the 
method of Lowe, except that the fraction (1 — P) of pairs which do not have their relative velocities stochastically 
reselected, are instead thermalized by a deterministic method. For each such pair, after steps (|5a|) . Ij5b() . a dissipative 
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force is calculated 

fij = few • fij) fij 

f?:=fF + f» 

Then, after steps ijSDl . (|5d|l . the momenta are corrected by 

Pl :=p t + At(l-f k /T)/ t D , 

where is a kinetic temperature estimated from relative velocities at the start of the timestep. Finally, the Lowe 
velocity reselection process is applied to the remaining fraction P of pairs as usual. In the present work, following 
Ref. [l8j , the parameter k is set to k = 0.3/At. As for the Lowe method, the key parameter is the stochastic 
randomization rate T which determines P. |l8j associate the deterministic part of their thermostat with the Nose- 
Hoover algorithm [2^, |2j| , but in practice it seems to be more closely related to that of Berendsen et al. (^3] , since 
the "friction coefficient" is not a dynamical variable but is instead directly proportional to the temperature factor 
(1 — Tt/T). This approach permits a high degree of tem perature control at low, or even zero, values of L, where the 
Lowe method is ineffective: indeed. Istovanov and Grootl report an order of magnitude improvement in 2V However, 
the links with a statistical ensemble, and with the DPD equations themselves, remain to be established. 

III. SIMULATION MODELS AND PARAMETERS 

By considering the equation of state and dynamical properties of water under ambient conditions, Groot and 
Warren Q established the standard parameter set that is commonly used in algorithm tests 0, 0, 0, : strength 
parameter a = 25, number density p = 3, stochastic impulse strength a — 3, and hence (for T = 1) 7 = 4.5. This 
is sometimes referred to as "Model B". This paper reports results for Model B, and also for variants with a = 1 and 
(7 = 6, using both the DPD and Peters methods. For the Stoyanov-Groot and Lowe algorithms, stochastic collision 
rates T = 0, 0.2, 4 and the maximum possible T = 1/ At (corresponding to P = 1, completely reselecting all velocities 
every timestep) are investigated. A system size of N = 250 particles in a cubic box with periodic boundary conditions 
is used throughout. 

For their tests of DPD simulations of membranes, Jakobsen et al. 26\ adopted the model of Shillcock and Lipowsky 
. Each lipid molecule has the form of a 7-bead chain HT6 in which repulsion parameters between hydrophilic "head" 
beads (H), hydrophobic "tail" beads (T), and "water" beads (W) are as follows: «ww = «hh = "tt = 25, auw = 35, 
«ht = 50, axw = 75. Successive lipid beads are connected by harmonic springs, and a simple bond-bending potential 
favours linearity: 

"stretcher,) = \ ^stretch (n j - £ ) 2 with fcstretch = 128, £q = 0.5 (6a) 

Ubend(rij,rjk) = fcbcnd(l - fij ■ fjk) with k hcnd = 20 . (6b) 

The present paper uses this model but follows Jakobsen et al. |2y| in choosing all stochastic impulse strengths to be 
the same, a — 3, in the DPD and Peters algorithms. For comparison, results obtained using the Stoyanov-Groot 
method with T = 0, and using the Lowe thermostat with T = 1/At, are reported. A smaller system is used here than 
in the previous papers 0,|2(| : Mipid = 100 lipid chains plus N watcl = 2500 water beads, making N — 3200 beads in 
all. A cuboidal box, with periodic boundaries, was employed. Transverse box dimensions were fixed so as to give an 
overall number density p = 3, and an area per lipid 2A/N\i V i<± = 0.617, close to the state of zero surface tension. The 
chains were pre-formed into a flat bilayer, and allowed to equilibrate before any measurements were made. 

The time unit r = r c ^Jm/k^T is defined in terms of the other fixed parameters. In these units, for the water 
model described above, Groot and Warren |3( recommend timesteps in the range At — 0.04. . .0.06, so as to achieve 
an accuracy of order 1% in Tt; the present paper investigates time steps in the range 0.005 < At < 0.06. Total 
simulation times for the water runs were t lun = 2000, and for the membrane runs were t run = 10000. Formulae for 
the Laplacians used to calculate the configurational temperature for this model are given in the Appendix. 

IV. RESULTS 

Figures H and show the measured temperatures for the four chosen algorithms, in simulating the water model at 
various choices of thermostatting strength. With an appropriate choice of parameters, it is indeed possible to achieve a 
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FIG. 1: Water simulations using (a) DPD algorithm with A = 0.65, (b) Peters scheme II. Temperatures as functions of timestep 
At. Solid lines: configurational temperature T c . Dashed lines: kinetic temperature Tk. Circles: a = 1. Squares: a — 3. 
Diamonds: a = 6. The lines are to guide the eye. Statistical errors are approximately the same size as the plotting symbols. 




FIG. 2: Water simulations using (a) Stoyanov-Groot thermostat, (b) Lowe thermostat. Temperatures as functions of timestep 
At. Solid lines: configurational temperature T c . Dashed lines: kinetic temperature T k . Downward-triangles: T = 0. Circles: 
r = 0.5. Squares: r = 4. Upward-triangles: V — 1/At (maximum possible). The lines are to guide the eye. Statistical errors 
are approximately the same size as the plotting symbols. 




kinetic temperature Tk within 1% of the desired value at a timestep At — 0.05; the Stoyanov-Groot method performs 
well in this regard at all values of T, and the Lowe thermostat is good at high values of T. The more recently-developed 
thermostats all seem less sensitive to variations in a or T than the basic DPD algorithm. However, it can be seen 
that the error in the configurational temperature T c is an order of magnitude worse than that in Tk, being around 
10% at At = 0.05, and this is not dramatically affected by the choice of algorithm: the "best" (Stoyanov-Groot at 
r = 0) reduces the error by 30-40% compared to the "worst" (all of the methods, at high stochastic thermalization 
rate). Also, for quite small stochastic damping T = 0.5 (i.e. P = 0.025 at At = 0.05), the Stoyanov-Groot method 
becomes similar in performance to the other methods. 
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FIG. 3: Membrane simulations using (a) DPD algorithm with A = 0.65, a = 3.0, (b) Peters scheme II with a — 3.0. 
Temperatures as functions of timestep At. Circles: H (head) and Te (tail) beads. Squares: Ti and T5 beads. Diamonds: T2, 
T3 and T4 beads. Triangles: W (water) beads. The lines are to guide the eye. Statistical errors are approximately the same 
size as the plotting symbols. 
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For the membrane simulations, temperatures for each of the bead types (species) are reported in Figs and 0] 
There are several points of interest. Firstly, the kinetic temperatures for different species are different, as already 
noted [2B|. The discrepancies between and T for the lipid beads are somewhat worse than in the case of simple 
water at the same timesteps, of order a few percent. The configurational temperatures are also different for the 
different species. For the lipid beads, they are all in much poorer agreement with the desired temperature, than for 
the water: at At = 0.05, T c is too high by 80-90% for some of the tail beads. This indicates a serious departure from 
equilibrium. Secondly, numbering the tail beads T\—Tq with Ti next to the head bead, the results for both Tk and T c 
divide naturally into groups: {H,Tg}, {T^Ts}, {T2,T3,T4}, and {W}. Within each group, the temperatures agree, 
almost to within the statistical error bars, and no attempt is made to distinguish them in the figures. This suggests 
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FIG. 5: Density and temperature profiles in the bilayer region for membrane simulations using the Peters thermostat. Top two 
panels show configurational and kinetic temperatures as functions of timestep. Solid lines: At — 0.05. Dashed lines: At — 0.03. 
Dot-dashed lines: At — 0.01. Bottom panel shows, for reference, the local density, measured at At — 0.01. Solid lines: H 
beads. Long dashes: T beads. Short dashes: W beads. 




z 



|26| that the intramolecular stretching and bending potential terms (which are the same within each group) are the 
dominant factors in the problems of maintaining equilibrium; in the following section, some support will be given to 
this view, but it is also argued that the intermolecular potentials make a contribution. Of the four algorithms studied, 
the DPD, Peters, and heavily-thermostatted Lowe method perform in a similar way to each other. The behaviour 
of the Stoyanov algorithm is unusual: although T c is quite close to T up to At = 0.03, 7k shows an instability. As 
At increases, the water temperature falls dramatically below T, while the other species become much too hot. This 
effect makes the simulation unphysical for At = 0.05. 

Density and temperature profiles are shown in Fig. [SI selecting the Peters thermostat as being representative. Jakob- 
sen et al. [2(| report a time step dependence of the density profiles, specifically a change in the (small) concentrations 
of head group beads in the interior of the bilayer. The present simulations do not confirm such a systematic effect; 
however, the transverse dimensions of the system are much smaller than those of Ref. [2(|, so the relevant statistics 
and timescales may prevent resolution of the effect here. The temperature profiles are as expected from Figs and 
0J with the most dramatic indicators of lack of equilibrium in the regions occupied by tail beads {T2,T3,T4}. A 
comparison of the different algorithms at fixed timestep At = 0.03 is given in Fig |SJ In the lipid tail region the 
Lowe thermostat (with the parameter T set to the maximum allowed value) gives the worst T c , about 20% too high, 
while the Stoyanov thermostat (with T = 0) gives the best, about 8% too high; however the performance order is 
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FIG. 6: Density and temperature profiles in the bilayer region for membrane simulations at At = 0.03 using the different 
thermostats. Top two panels show configurational and kinetic temperatures. Solid lines: Lowe thermostat with F = I/At. 
Long dashes: Peters thermostat with a = 3. Short dashes: DPD algorithm with a — 3. Dash-dotted line: Stoyanov-Groot 
thermostat with V = 0. Bottom panel shows, for reference, the local density, measured at At — 0.01. Solid lines: H beads. 
Long dashes: T beads. Short dashes: W beads. 




reversed for with the Lowe thermostat (by construction) being perfectly accurate and the Stoyanov thermostat 
being overdamped by about 5%. The behaviour in the water region does not give a good indicator of the lack of 
equilibrium within the membrane. 

It is of interest to compare these profiles with the pressure tensor results which were also prev iously found to 
exhibit an anomaly pf3 | . The pressure profiles are calculated in the Irving-Kirkwood convention |3l| using a method 
essentially the same as that described by Goetz and Lipowsky 32] . The results for the Peters thermostat are shown 
in Fig. [7| The component Pn normal to the bilayer is not constant, contrary to what must hold for a system at 
equilibrium |26| . For At = 0.05 the error is comparable in magnitude to the physically interesting quantity, the 
difference Pn ~ Pr between normal and transverse components, whose integral gives the surface tension. There is 
some cancellation of errors, since the anomaly appears in both components 01 however, such an error is clearly 
undesirable. The point to be made here is that the configurational temperature profiles of Fig. [3] are echoed in the 
pressure profiles of Fig. [7| 
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FIG. 7: Density and pressure profiles in the bilayer region for membrane simulations using the Peters thermostat. Top panel: 
normal component Pn as a function of timestep: Solid lines: At = 0.05. Dashed lines: At = 0.03. Dot-dashed lines: At = 0.01. 
Middle panel: surface tension integrand, Pn — Pr, measured at At = 0.01. Bottom panel: local density, measured at At = 0.01. 
Solid lines: H beads. Long dashes: T beads. Short dashes: W beads. 




-10 -5 5 10 

Z 



V. ANALYSIS AND DISCUSSION 

It is evident that the configurational temperature gives a different perspective, from the kinetic temperature, on 
simulation equilibrium, even for simple fluids. Also the intramolecular bond stretching and bending potentials used in 
membrane simulations have a direct effect on the stability of the algorithm, which can be estimated via T c , amongst 
other indicators. 

Some insight into the information contained in T c may be obtained by considering a very simple model: the 1-D 
harmonic oscillator, with hamiltonian 

U = p 2 /2m + \kx 2 = K(p) + U(x) . 

It is known [s^l that the velocity Verlet algorithm conserves exactly a "shadow hamiltonian" which we may write 
in the form 

H x {x,p) = p 2 /2m+\kx 2 {l -<f>)= K{p) + U*(x) 
where <f> = j(k/m)At 2 . 
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TABLE I: Laplacian Vf U for beads in lipid chain calculated for isolated chain at T = 0, as an equilibrium simulation average 
at T — 1 for an isolated chain, from the full condensed phase simulation at T = 1, and from a fit to the measured T c , Fig. 0J:>, 
assuming a locally isotropic harmonic potential. Figures in parentheses represent the estimated error in the last digit. 





single 


chain 


condensed phase 


Bead 


T = 


T=l 


T = 1 


Fit 


H 


288 


273(2) 


491(1) 


1157(1) 


Ti 


1056 


853(6) 


1132(1) 


1893(1) 


T 2 


1216 


979(5) 


1272(1) 


2272(1) 


T 3 


1216 


986(7) 


1266(2) 


2327(1) 


T 4 


1216 


972(8) 


1264(2) 


2277(1) 


T 5 


1056 


865(8) 


1116(2) 


1884(1) 


T 6 


288 


267(3) 


438(1) 


1123(1) 


W 






122.8(1) 


479(1) 



A dynamical algorithm which alternates periods of velocity Verlet with velocity reselection from the canonical ensemble 
at temperature T (analogous to the "Lowe" algorithm of DPD) will actually generate momenta from the correct 
Maxwell-Boltzmann distribution exp(— K./k B T) and configurations from the "shadow distribution" exp(— U^/k&T). 
An immediate consequence is 

\k B T = <p 2 /2m) = (\kx 2 (l -0)) . 
Since the force / = — kx and d 2 7i/ dx 2 = k, the kinetic and configurational temperatures are 

k B T k = (p 2 /m) = k B T (7a) 

, rp (P) k B T 

kBTc = (d 2 n/dx 2 ) = w • (7b) 

For a 3D isotropic harmonic oscillator U = ^k\r\ 2 , f = —kr the same formula applies, with k = ^W 2 U, and it is 
straightforward to extend it to the anisotropic case. The effect of an over-large time step is to reduce the effective 
force constant dictating the configurational distribution, increase the mean squared displacement, and (in proportion) 
increase the mean squared force which gives the value of T c . The above equation is a reasonable first approximation 
to the curves shown in Figs ^ [5] For the lipids, assuming that the intramolecular bond potentials are approximately 
harmonic, a given fractional error in T c for each bead may approximately translate into a corresponding fractional 
error in mean-squared displacement which, when compounded along the chain, becomes a more significant error in 
overall flexibility. 

It is tempting to try to estimate these contributions to T c , for the isolated lipid chain. In the harmonic approxima- 
tion, the appropriate force constants may be evaluated, taking the lipids to be in their minimum energy configuration, 
and counting all the bonds in which a bead participates. In Table ^ the sum (over x, y and z directions) of these 
predicted force constants (for fc s t re tch = 128, £q — 0.5, fcbcnd = 20) are compared with measured average Laplacians 
from simulations of a single chain at T = 1, and simulations of the full fluid system at T = 1, both conducted at 
short timestcps. The figures partly explain the close similarity in configurational temperatures within the groups of 
beads mentioned in the previous section, in terms of intramolecular potentials. There is some reduction in the average 
Laplacians due to fluctuations of the single chain at T = 1 compared with T = 0. However, in the full condensed 
phase simulation, this is more than compensated by an intermolecular contribution of order 200-250 in these reduced 
units, so part of this coincidence is due to the similar environments, and identical repulsion parameters between like 
beads, for this particular model. 

The simulated fluids are not actually harmonic, and so the above analysis should not be over-interpreted. The final 
column of Table [i] shows the results of fitting eqn (7E)| to the Lowe thermostat results shown in Fig.^J). The fits are 
excellent, but the fitted parameters are all significantly larger than would have been expected. Partly, this can be 
explained by the anisotropy of the local harmonic potential: the stiffest potentials will dominate the breakdown of the 
algorithm. However one should also remember that the real shadow hamiltonian is more complicated than a simple 
harmonic potential, and will not split in the same way into a kinetic part plus a renormalized configurational part. 

Some effort has been made to develop practical methods to estimate the shadow Hamiltonian for complex systems 
|35l l36j| . Knowledge of this function may help not only the monitoring of algorithm performance, but also the 
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implementation of more efficient algorithms such as hybrid Monte Carlo |3j|. This method differs from DPD in the 
rejection of moves so as to guarantee that the desired ensemble is sampled. Returning to the points made by Hafskjold 
et al. |27| . it is notable that smoothness of the potential is a highly desirable feature in the accurate construction of 
shadow hamiltonians. To illustrate this, Engle et al. |36j| considered a simple harmonic oscillator potential, modified by 
splitting into two halves at the minimum, with a flat piece inserted between them. The consequent lack of smoothness 
was shown to degrade conservation of the shadow hamiltonian. Such a potential is a possible one-dimensional model 
system for the numerical aspects of DPD dynamics: the oscillation back and forth represents the motion of a DPD 
bead out of the repulsive quadratic potential of a neighbour, and into the repulsive range of a different neighbour. 
In the absence of an improved algorithm to handle the DPD conservative forces, using a shorter timestep for them 
seems the best approach |l7j . 

VI. CONCLUSIONS 

The essential conclusion of this paper is that stiffer potentials, and stronger conservative forces, re quir e shorter 
timesteps. This is hardly a new point, and is increasingly being recognised in the DPD community |l7l Efj l27j . 
where the principal danger is that the effects of an over- long timestep may be disguised by the thermostat. It is 
recommended here that the configurational temperature T c be measured and reported in DPD simulations, along 
with other indicators that the system is not at equilibrium, so as to address this danger. T c is easy to calculate, and 
appears to be sensitive to some of the artifacts that appear in other properties such as the pressure tensor profile in 
membranes. 

The degree of accuracy of sampling required will, no doubt, vary from one application to another. For the standard 
water model, Groot and Warren recommended timesteps At rj 0.05, so as to achieve an accuracy of order 1% in 
Tjj. Applying the same criterion to T c would lead to a revised recommendation, At rj 0.015. Modern thermostats 
are dramatically better at controlling Tk, and give modest improvements in T c . Stoyanov and Groot 0] note, for 
instance, that their new method keeps Tk within the desired 1% range even for At = 0.09; however, such a timestep 
would result in an error in T c between 24% and 44% depending on T. Also, for multicomponent systems such as lipid 
membranes, there may be substantial imbalances in the distribution of Tk and T c between the different species. Such 
large discrepancies indicate lack of equilibrium, i.e. configurations are being sampled incorrectly. 

It would, no doubt, be possible to devise thermostats that control the value of T c as well as Tk, and that act 
separately on the different types of bead. In such a case, it becomes important to measure other, independent, 
quantities as a function of time step, and to bear in mind that no single quantity acts as an indicator of equilibrium. 
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Appendix 

The calculation of the Laplacian for each of the potential terms used in this study proceeds as follows. The DPD 
repulsive potential between i and j, with m = Ti — rj, takes the form (assuming r.y < r c ) 




The gradients (the negatives of which are the conservative forces of eqn (|3a(l ) are 

Viit = — VjU 

and the Laplacians are 

V^m = Vftt 

The formulae for the bond stretching potential of eqn l|tja[) are identical, with the replacements r c — > £ and a 

^strctch^O- 
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For the angle-bending potential of eqn i|6b[l between successive beads i, j, k, define cy = r,j ■ m, Cj k = rjk ■ rjk, 
"ijk = r ij ' r jk- The gradients of these scalar products with respect to individual positions are 



r 3k 



^ iCij — ^^ij 


ViCjfe = 








^ jCjk — 2Vj k 


^ 'jCijk 


r «i — 


VkCi-j = 


VjfeCjfc = -2r,- fe 




~ r ij 


are 








V?Cy = 6 


V?c jft = 




= 


Vf Cy = 6 


V|cjfc = 6 




= -6 


V 2 feCy = 


VfeC jfc = 6 




= . 



The potential has the form Mbend = &bend(l — C), where 



C = CQs6 ijk = c tJ 1/2 c jk /2 c ijk 



The necessary gradients are 



-1/2 -1/2 / Cijfe 

ViC = (v ) ^)^ ] ^ 1/2 (v J ^% t + c , I V /2 (v )Cl]l ) 
= c -v 2c -v 2 ( (l + ^_ (l + ^ 

VfeC = (V fc c 4 r. 1/2 )cT fc 1/2 C y7c + Cy 1/2 (V fe c~ fc 1/2 )c 4ife + Cy 1/2 c7 fe 1/2 (VfcC 4ifc ) 

- r- 1 ^ -1/2 f Cijk 

- ij > k \c jk Tjk 
The Laplacians are given by 

+ 2c^ /2 (V lC j k 1/2 ■ V iCijk ) + 2cJ k 1/2 (V t c l3k ■ V 4 c-. 1/2 ) + 2c l]k (V t c^ /2 ■ V.c" 1 ^ 

3/2^-1/2^ 

U 



2c„„- ' Cjj, ' c^k 



V 2 C = (Vfo^cJ^ + c^^^ 

+ 2c^ /2 (V jC ; k 1/2 ■ V jCijk ) + 2cJ k 1/2 (V jCl]k • V jCi r. 1/2 ) + 2c ijk (V j cT j 1/2 ■ V 3 cJ k 1/2 ) 

= ~^ c ij ^ c jk ^ { C ijk + Cijk(Cij + Cjk) + CijCjk) 
V 2 r _ / V 2 -1/2N -1/2 -1/2/WJ -1/2N -1/2 -1/2^2 N 

+ 2c^ /2 (V k cj h 1/2 ■ V k c ijk ) + 2cj k 1/2 (V k c l]k ■ V k c~ 1/2 ) + 2c ljk (V k c^ /2 ■ V k cj k } /2 ) 



.-1/2^-3/2^ 



2c^ A ^jk Cijk 



All of the above results should be multiplied by (— fcbcnd) to give the gradients and Laplacians for the bending potential 
of eqn 

It is worth emphasizing that the dissipative and random terms in the DPD equations (J2J) play no role in the 
calculation of Laplacians and squared forces which go into the configurational temperature, eqn 1(1)1, 
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